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FLOW OF ELECTRONS THROUGH A RAREFIED GAS 
BETWEEN PLANE PARALLEL ELECTRODES 
by Peter M. Sockol 
Lewis Research Center 

SUMMARY 

The integral equation approach to the solution of the Boltzmann equation has been ap- 
plied to the flow of electrons through a rarefied gas between plane parallel electrodes. A 
perfect Lorentz gas model is assumed. Self- consistent-field solutions of the Boltzmann 
equation with Poisson’s equation have been obtained for ratios of electrode spacing to 
mean free path of 0, 0. 5, 1, and 2. Density and potential curves together with current- 
voltage characteristics are presented. 


INTRODUCTION 

The flow of charged particles between an electrode and an ionized gas has attracted 
the attention of several investigators since the early work of Langmuir. Most of these 
studies have been concerned with the theory of electrostatic probes (refs. 1 to 9). Some 
work, however, has been done in connection with the analysis of thermionic converters 
(refs. 10 and 11). All these treatments have involved one of two approaches. In the first, 
either collisionless or diffusion equations (depending on the density range of interest) 
have been applied continuously from the electrode to the undisturbed portion of the gas or 
to another electrode. In the second, the space has been divided into two or more distinct 
regions; appropriate equations, either collisionless or diffusion, have been applied in 
each region, and more or less arbitrary conditions have been used to match solutions at 
the boundaries between regions. 

In situations where collisions are important, both of these approaches are question- 
able. There will inevitably be a transition region, however brief in extent, in which 
neither collisions nor particle dynamics may be considered as dominant. From the ex- 
periences of investigators studying this regime in ordinary gas dynamics, it would seem 
that no simple solvable model is generally valid here. Hence, in the present work, an 



attack is made on the problem by means of the Boltzmann equation. Rather than treating 
the problem in all its complexity, however, a simpler problem was treated as completely 
as possible, namely, the flow of electrons through a rarefied gas between parallel elec- 
trodes. Taking one of the electrodes to be thermionically emitting and the other cold al- 
lows some of the aspects of flow between an ionized gas and a cold electrode to be re- 
tained. The next steps would be the inclusion of positive ions and then the substitution of 
an ionized gas for the hot electrode. 

Possible methods for solving the Boltzmann equation for such a boundary value prob- 
lem are provided by the half-range polynomial method of Gross, Jackson, and Ziering 
(ref. 12) and the discrete ordinate method of Chandrasekhar (ref. 13). The inclusion of 
a spatially varying electric field, however, makes these methods intractible. Instead, 
the original integro- differential equation is converted to a pure integral equation and is 
solved numerically. The method is similar to that used by Willis to solve the Krook 
kinetic model for shear flow (ref. 14). 

The particular problem considered herein has been treated recently by Goldstein, 
using Monte Carlo techniques (ref. 15). While the present method is perhaps somewhat 
more difficult to apply, it offers the advantage of obtaining information about the distribu- 
tion function itself. 


FORMULATION 

The problem under consideration is that of the steady flow of electrons between infi- 
nite plane-parallel electrodes. The Boltzmann equation for the electron distribution 
function f may be written as (ref. 16) 


vcosel + iLel + ^i — 

3x x \ 3v v 3 cos 6 > 



( 1 ) 


where x, v cos d, and a. x are, respectively, distance, velocity, and acceleration nor- 
mal to the electrodes. For electrons interacting with the space charge field 


m dx 


( 2 ) 


where V is the electrostatic potential. Bars have been used to designate dimensional 
quantities. 

For electrode separations that are not too large compared with the mean free path, 
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the loss of energy by electrons due to elastic collisions with atoms is negligible. In addi- 
tion, all inelastic and electron- electron collisions are neglected. The collision term in 
equation (1) thus reduces to that for a perfect Lorentz gas (ref. 17): 

(st") = Nv / [f(x,v, 0’) -7(x,v, e)]cr(v,x)do>» (3) 

c 

where N is the atom density, cr the differential scattering cross section, x the scatter- 
ing angle, and do)’ the element of solid angle. Primes refer to values after a collision. 
The number density and the particle current are given in terms of f by 


and 


— r°° S ' 71 — 

n = 2 tt / / f sin 9 

Jo Jo 


d 9 v 2 dv 


J = 277 


r °° 
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—3 — 
f cos 9 sin 9 d0 v dv 
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( 5 ) 


Equations (1) to (4) are to be solved self- consistently with Poisson's equation 



= 477 en 


( 6 ) 


The boundary conditions on the distribution function f are of a rather peculiar na- 
ture; f is not completely known on one surface in phase space, such as that specified 
by x = 0. The number of particles of a given velocity returned to a surface by back 
scattering from the gas is not known a priori. Instead, the emergent half of f is known 
at each of two bounding surfaces, x = 0 and x = L. Hence, with half- Maxwellian emis- 
sion of electrons at x = 0 and complete absorption of electrons at both x = 0 and x = L, 
the boundary conditions on f and v are 


x = 0, cos 9 > 0, f = 2n [• — 2) — 

°\27TkT 


\3/2 -mv 2 /2kT 


x = L, cos 9 < 0, f = 0 


V(0) = 0, V(L) = V 


( 7 ) 
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where n Q and T Q are the density and temperature of the emitted electrons, L the elec- 
trode spacing, and V c the potential of the cold collector. 

Before proceeding further, it is necessary to assume that the differential scattering 
cross section remains finite for small scattering angles x* Then the collision integral, 
equation (3), may be separated into two finite integrals. The cross section can then be 
approximated by a finite sum of Legendre polynomials in cos x> and then in turn, ex- 
pressed in terms of cos 9 and cos 9 r by means of the addition theorem for spherical 
harmonics. This procedure would assume more information than is generally available 
on electron-atom cross sections. Hence, in the present work, a is assumed independent 
of x- Equation (3) is rewritten as 


f—\ = -X ^vf+ — X - ^v f f(x, v, 9 1 ) sin 9 1 d9' 
\6t/ 2 SO 


( 8 ) 


where the mean free path X is given by 


X - ^ = 47tNct 


(9) 


and the element of solid angle do?' has been expressed as 


da)' = 2n sin 9' d0* 


( 10 ) 


Dimensionless variables are now introduced: 


x — m 
x =— , v = v 

L \2kT 


a/2 


, p = cos 9, <p = 


eV 

kT 
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Equations (1), (2), and (8) are combined as 
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(11) 


3f 
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and equations (4) to (7) become 



<^=Cn 

j. .2 


- 2 1 

x = 0, p>0, f=e v 


x=l, p<0, f = 0 > 


<p(0) = 0, <p( 1) = (p c 


( 12 ) 


(13) 


(14) 


(15) 


t 


The transformation of equation (11) to an integral equation is accomplished by means 
of a change of variables from (x, v, ju) to (x, w, u) where 


w = v - cp(x) 
u = - (p(x) 


(16) 


The nondimens ional total energy of an electron is denoted by w, and u is the nondimen- 
sional portion of the total energy associated with the motion in the x- direction; thus, w 
and u are constants of the motion between collisions. Since f is a double valued function 
of u, f is denoted by f* for p. ^ 0. Finally, with note taken of the boundary conditions 
(15), g* is defined by 


f ± = e" 


w g± 


(17) 


Then equations (11) to (13) become 
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( 18 ) 


= TKCl'(x)g ± ± - 
9x 2 


K S(x,w) 
Vu + <p(x) 


n = A f e _w S(x,w)dw 


(19) 


with 


/*°° /*W 

J = / e" w / (g + - g")du dw 

y-<p( X ) ^-<?(x) 


or(x) 


¥ u + 


<p(x) 


cp(x) 


( 20 ) 


( 21 ) 


/*W 

S(x, w ) = - / (g + + g~) dU 

J-<p(x) Vu + (0(X) 


(22) 


The boundary conditions (15) become 


x = 0, g + = 1 




x = 1, g" = 0 > 


<P( 0) = 0, cp(l) = <p 


(23) 


Equation (18) is formally integrated to give 


± A ± 
l - A exp 


f x 

i f x 

f x 

T / /ca(z)dz 

±- / S(y,w)exp 

T / KQ!(z)dz 

5- 

0 H- 

1 

2 X ± Vu + <p(y) 
0 

U J 


dy 


(24) 


The first term in equation (24) is the contribution to the distribution function from parti- 
cles that have not suffered collisions; it decays exponentially as the ratio of the path 
length to the mean free path. The second term represents the contribution from parti- 
cles scattered into a particular ’’volume" element in phase space. The A ± and x* 
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are arbitrary functions of w and u. They are 
to be determined from the boundary conditions 
(23) together with continuity considerations at 
turning points for electrons in the electric field 
(see following section). Note that, in general, 
k will be a function of (w + q> (x)) . 

Equation (24) is an integral equation for 
g*, as S(x, w) is functionally dependent on g. 

A single linear integral equation for S(x, w) may 
be obtained by eliminating g ± between equations 
Figure 1. - Diagram of u,x plane. (22) and (24). 



ELECTRON FLOW WITH POTENTIAL MINIMUM 


The equations of the previous section are now applied to a situation in which a poten- 
tial minimum occurs within the interelectrode space. All other real situations can be 
obtained as special cases of this one by defining the potential minimum to be the point of 
lowest potential in the system even when this occurs at one of the electrodes. 

Consider now the diagram (fig. 1) of the u, x plane for a given value of w. From 
equation (16), it is seen that u > -cp(x), and thus the surface u = -cp(x) is a lower bound- 
ary for the phase space. At this surface, electrons in region I are reflected and transfer 
from g + to g"; the opposite is true for region n. Conservation of electrons requires 
that g + = g~ on this surface. The complete conditions on g ± are as follows: 


For 0 < x < x^, -cp(x) < w < -<P m > ~<P( X ) — u < w (region I), 

x = 0, g + = 1 > 

► 

X = x t , g- = g-^ 

For 0 < x < x m , -<p m < w < oo, -cp(x) < u < -<p m (region I), 

x = 0, g + = 1 > 

> 

X = X, , g" = g + 

L J 


( 25 ) 


(26) 
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For x m < x < 1, -cp m < w < -<p(x) < u < -cp m (region II), 


x = x t , g + = g' 
x = 1, g" = 0 


For 0 < x < 1, -<p m < w < < u £ w (region m), 

x = 0, g + = 1 > 

► 

x = 1, g" = 0 > 


where 


w + = 0 
u + <p(x t ) = 0 


( 27 ) 


(28) 


(29) 


The points Xj.(u) are turning points for electrons of ,f energy” u in the electric field. 

Condition (25) is the same as condition (26) and is not considered separately. Appli- 
cation of the conditions to equation (24) gives the following equations: 


For 0 < x < x , -cp(x) < u < -cp m (region 1), 


g + = exp 


/ x nX 

kch(z) dz + - / K - 

2 Jo Vu + <p(y) 


S(y,w)exp 


f 


fca(z)dz 


dy 


g = exp 


f x t /* x t 

- I Ka(z)dz - / /co(z)dz 

_ JO Jx 

.i n K 

2 Jo 

1 n x 

2 J x Vu + <p(y) 


V (30) 


S(y,w)exp 


S(y,w)exp 


r x t /* x t 

I /caf(z)dz - / 

t/y «/x 

r 

I /ca(z)d 

Jx 


dy 


J 
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For x m < x < 1, -cp(x) < u < -cp m (region II), 


.... 

2 J x Vu + <p(y) 


S(y,w)exp 


r y r 

■ I ko! ( z)dz - I KO!(z)dZ 

e/X. t/X, 


dy 


+i r_£_ 

2 X Vu + «p(y) 


S(y,w)exp 


[/ 


k ot (z)dz 


dy 


y (3D 


f- 


S(y,w)exp 


+ <p( y) 


/ 


k(x(z) dz 


dy 


For 0 < x < 1, -<p m s u ^ w (region III), 


g + = exp 


f x 

1 f x 

f* 

- 1 K 0 l(z)dz 

+ - / S(y,w)exp 

- / k ot (z)dz 

[Jo J 

2 Jq Vu + <p{ y) 

Jy J 




dy 


g 


-i r - * 

2 X Vu + <p(y) 


f (32) 


S(y,w)exp 


/ 


KQ!(z)dZ 


dy 




Examination of equations (30) to (32) shows that there are discontinuities in g" be- 
tween regions I and m and in g + between regions n and m. These discontinuities are 
permissible, however, as u = -<P m is a characteristic surface for the original integro- 
differential equation (11). 

Substitution of equations (30) to (32) into equation (22) , which defines S(x, w) , and 
interchanging the order of the u and y integrations gives a Fredholm integral equation 
for S(x,w). Care must be taken in interchanging the order of integration, as x^. is a 
function of u. For example, 




du dy 
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where 


u + (p(x t ) = 0, b + <p(x b ) = 0 


The integral equation may be written as 


S(x, w) = F(x, w) + - / kH(x, y, w)S(y, w)dy 

2 */0 


where 


x w , -<p(x) < w< 


X - w 


if 

c/a x 


F(x, w) = — / exp - I /co^dz 


/u + <p(x) 


+ — f exp- f koi(z) 6 z - /* KOf(z)dz — uu (35) 

2 X, X X 


■if 


H(x, y, w) = — / exp- / «Q!(z)dz 


/u + <p(x) \u + «p(y) 


/ b 

expj-j^y K«(z)dz - 


l r 

| I KOI 
I t/x. 


(z)dz u u (36) 

Vu + ^ (x) ^ (y) 


w, -cp(x) < w < -<p 


.-«W _< ^m — w °° 
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(38) 


r -(p(x), 0 < x < x. 


a f = 


m 


-cp . x < x < 1 
^ ^m’ m 


\-(p(x), cp(x) < cp( y) 

a h = ^ -<p(y)» <p(y) ^ <p( x ) 


~b, y — x m — x or x - x m - y 


(39) 


From equations (36) and (39), it may be seen that H(x, y,w) is symmetric in x and y. 

For x > x m , the distribution function vanishes for w < -<P m - Hence, equation (19) 
for the density may be rewritten as 


— f e” w S(x,w)dw (40) 


The current is most conveniently evaluated at the emitter (x = 0). Only those elec- 
trons whose energy w > -<P m can contribute to the current; hence, substitution of equa- 
tions (30) and (32) for g* into equation (20) and interchanging the order of the u and y 
integrations give 


J = e 


^“(1 - <p m ) - f e” w f~ (Pm exJ-2 f Xt K a(z ) dz 

J-<p m J ° L y ° 


du dw 


- f e“ w f /cS(y, w)F(y, w)dy dw (41) 

J-q> 

Y m 


NUMERICAL INTEGRATION 

The solution of equation (33) presents certain difficulties as the kernel H(x, y, w) is 
singular for y = x. The following method, however, proved to be quite satisfactory. 

On deleting the w, equation (33) may be written as 

S(x) = F(x) + — f 1 /c(y)H(x, y)S(y)dy (42) 

2 JO 
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(43) 


Now as y — x. for x + x„ 

m 

H(x, y) - - log(jx - y |) + const + • • • 

This fact is utilized to replace equation (42) by a system of linear algebraic equations. 
Equation (42) is rewritten as 

S(x) = F(x) + - L(x)S(x) +- f 1 /c(y)[log(|x - yhs(x) + H(x, y)S(y)ldy (44) 

2 2 J Q 


with 


L(x) = - f 1 /c(y)log(|x - y|)dy 

*/o 


(45) 


The integral in equation (44), whose integrand is everywhere finite, is replaced by a 
quadrature sum 



n 

f(x)dx (a Yj a j f ( x j) 

j=l 


(46 


and the equation is evaluated at each of the quadrature points to give the desired set of 
algebraic equations 



(i = 1, 


> n) 


(47) 


where 


Si = S(x.) 
Fj = F(x.) 
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« *i) 


a j K ( x j) H ( x i > x j) 


L(x.) + a i K(x i )H £ (x i ) + Yj a k K ^ log (l x i " x kO 0 = i) 

k=l 

k*i 


(48) 


and 


H £ (x) = lim[log (| x - y | ) + H(x, y)] 
y-x 


(49) 


This limit can be found numerically. 

The set of equations (47) can be solved by standard techniques. As H(x, y) is sym- 
metric in x and y, the set may be rearranged into a form with a symmetric matrix. 
Care must be taken in the choice of the quadrature formula (46). The function F(x) has 
a large negative derivative (in some cases infinite) at x = 0. To achieve sufficient accu- 
racy from the approximation (eq. (46)), the quadrature points must be grouped closely 
together near x = 0. In the present calculations, this was accomplished by subdividing 
the range of integration into two parts and applying Gaussian quadrature formulas in each 
section. The point of subdivision was taken quite close to x = 0. 

The complete calculational procedure can be outlined as follows. An initial potential 
distribution was assumed. Equations (47) were solved for S(x,w) for selected values of 
w, and the density n(x) was obtained from equation (40) . By means of various curve- 
fitting techniques, a Chebyshev series was fitted to n(x), and the potential cp(x) was ob- 
tained from Poisson's equation (14) by analytically integrating the Chebyshev series. 

This Picard iteration sequence was continued with the new potential for cp(x) until con- 
vergence was achieved. The current was then obtained from equation (41) with values 
from the final iteration. Solutions from the collisionless case were used as initial iter- 
ates for cases with low values of k. As collisional results became available, they were 
used as initial iterates for cases with higher values of k. 


RESULTS 

The equations have been solved on a high-speed computer for the case of rigid 
sphere scattering ( k constant). The particular solutions were for k = 0, 0. 5, 1, and 2 
with C (eq. (14)) = 50. An emitter temperature of 1900° K, a spacing of 10 microns. 
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0 .2 .4 .6 .8 1.0 


Dimensionless distance, x/L 

Figure 2. - Variation of density and potential in interelectrode space. 

Initial slope, p'(o) = 0; space charge parameter, C = 50 (eq. (14)). 

and an emission current of 10 amperes per square centimeter would give a value for C 
of about 50. 

Figure 2 shows the variations of density and potential, respectively, with x and k 
for the potential with an initial slope cp'( 0), of 0. The steep initial slopes in the den- 
sity curves reflect the similar behavior in the function F(x,w) of equation (35). 

Figure 3 shows the current- voltage characteristics obtained as the result of many 
solutions of the combined equations. The agreement between the present results and 
those of Goldstein (ref. 15) is of a peculiar form. While the points appear to lie on the 
same curve, attention must be given to the manner in which these points were obtained. 
In both calculations, cp and cp* were specified at x = 0; nevertheless, Goldstein's 
points lie somewhat to the left of those obtained in this work for the same conditions. 
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Dimensionless current, J 



Figure 3. - Effect of electron scattering on current-voltage characteristics. Space charge parameter, C - 50 (eq. (14)). 


This discrepancy is well outside the standard deviation of the Monte Carlo calculation at 
the low end of the current -voltage characteristic. The difference is probably due to 
numerical difficulties in one or both of the treatments. 

CONCLUSIONS 

The integral equation approach to the solution of the Boltzmann equation has been 
demonstrated for a problem of electron flow through a rarefied gas over distances com- 
parable to the electron mean free path. The method should be capable of extension to 
more complicated problems including the flow of charged particles between an electrode 
and an ionized gas. In addition, as the method provides information about the distribution 
function itself, it could be used as the basis of an investigation into approximate methods 
of solving these types of problems. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 7, 1966. 
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APPENDIX - SYMBOLS 


cljj 

integration limit (eq. (38)) 

V 

velocity 

a h 

integration limit (eq. (39)) 

w 

dimensionless energy variable 

a i 

quadrature weight (eq. (46)) 


(eq. (16)) 

J 

3* 

acceleration 

X 

distance 

b 

integration limit (eq. (37)) 

x j 

quadrature point 

C 

space charge parameter 

x t’ x w 

turning points (eq. (29)) 


(eq. (14)) 

X 1 

integration limit (eq. (42)) 

e 

electronic charge 

a(x) 

function (eq. (21)) 

F(x,w) 

inhomogeneous function if 

e 

polar angle in velocity space 


integral equation 

K 

ratio of spacing to mean free 

f(x, v,0) 

distribution function 


path 

G ij 

matrix (eq. (48)) 

X 

mean free path 

g(x,w,u) 

transformed distribution 

M 

cos 6 


function 

ct(v,x) 

collision cross section 

H £ (x) 

limit function (eq. (49)) 

(p(x) 

dimensionless potential 

H(x,y,w) 

kernel of integral equation 

X 

scattering angle 

J 

dimensionless current 

CO 

solid angle 

L 

electrode spacing 

Subscripts: 


L(x) 

logarithmic integral (eq. (45)) 

c 

at collector 

m 

electronic mass 

m 

at potential minimum 

N 

atom density 

Superscripts: 

n(x) 

electron density 

+ 

flow toward collector 

% 

electron emission density 

- 

flow toward emitter 


S(x, w) source function (eq. (22)) 

T q electron emission tempera- 

ture 

u dimensionless energy varia- 

ble (eq. (16)) 

V(x) potential 
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